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Abstract. This paper develops and analyzes an interior penalty discontinuous Galerkin (IPDG) 
method using piecewise linear polynomials for the indefinite time harmonic Maxwell equations with 
the impedance boundary condition in the three dimensional space. The main novelties of the proposed 
IPDG method include the following: first, the method penalizes not only the jumps of the tangential 
component of the electric field across the element faces but also the jumps of the tangential component 
of its vorticity field; second, the penalty parameters are taken as complex numbers of negative 
imaginary parts. For the differential problem, we prove that the sesquilinear form associated with 
the Maxwell problem satisfies a generalized weak stability (i.e., inf-sup condition) for star-shaped 
domains. Such a generalized weak stability readily infers wave-number explicit a priori estimates 
for the solution of the Maxwell problem, which plays an important role in the error analysis for the 
IPDG method. For the proposed IPDG method, we show that the discrete sesquilinear form satisfies 
a coercivity for all positive mesh size h and wave number k and for general domains including non- 
star-shaped ones. In turn, the coercivity easily yields the well-posedness and stability estimates (i.e., 
a priori estimates) for the discrete problem without imposing any mesh constraint. Based on these 
discrete stability estimates, by adapting a nonstandard error estimate technique of I10| . we derive 
both the energy-norm and the L-^-norm error estimates for the IPDG method in all mesh parameter 
regimes including pre-asymptotic regime (i.e., k^h > 1). Numerical experiments are also presented 
to gauge the theoretical results and to numerically examine the pollution effect (with respect to k) 
in the error bounds. 
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1. Introduction. This paper develops and analyzes interior penalty discontin- 
uous Galerkin (IPDG) methods for the following time harmonic Maxwell problem: 



where ft cz R'^ is a bounded domain with Lipschitz continuous boundary dfl and of 
diameter R. v denotes the unit outward normal to dil, i := V^lj the imaginary unit, 
and Et = (f x E) x i/, the tangential component of the electric field E. k, called wave 
number^ is a positive constant and A > is known as the impedance constant. (jl.2p is 
the standard impedance boundary condition. Assume that g • f = 0, hence, gT = g- 
Problem (|l.ll) - (jl.2[) is a prototypical problem in electromagnetic scattering (cf. [6] 
and the references therein) and has been used extensively as a model (and benchmark) 
problem to develop various numerical discretization methods including finite element 
methods [T71[21] and discontinuous Galerkin methods [TH (TSl UHl [3 [H] , and to develop 
fast solvers (cf. [52] and the references therein). The above Maxwell problem with 
large wave number k is numerically difficult to solve mainly because of the following 
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two reasons. First, the large wave number k implies the small wave length £ := 27r/fc, 
that is, the wave is a short wave and very oscillatory. It is well known that, in 
every coordinate direction, one must put some minimal number of grid points in each 
wave length in order to resolve the wave. Using such a fine mesh evidently results 
in a huge algebraic problem to solve regardless what discretization method is used. 
Practically, "the rule of thumb" is to use 6 — 10 grid points per wave length, which 
means that the mesh size h must satisfy the constraint hk < 1. To the best of our 
knowledge, no numerical method in the literature has been proved to be uniquely 
solvable and to have an error bound under the mesh constraint hk < I for the above 
Maxwell problem. Moreover, numerical experiments have shown that under the mesh 
condition hk < 1 the errors of all existing numerical methods grow as the wave number 
k increases. This means that the error is not completely controlled by the product 
hk and it provides strong evidences of the existence of so-called "pollution" in the 
error bounds. It is known now [5] that the existence of pollution is related to the 
loss of stability of numerical methods with large wave numbers for the scalar wave 
equation, which is also expected to be the case for the vector wave equations. Second, 
for large wave number fc, the Maxwell operator is strongly indefinite. Such a strong 
indefiniteness certainly passes onto any discretization of the Maxwell problem. In 
other words, the stiffness matrix of the discrete problem is not only very large but 
also strongly indefinite. Solving such a large, strongly indefinite, and ill-conditioned 
algebraic problem is proved to be very challenging and all the well-known iterative 
methods were proved numerically to be either ineffective or divergent for indefinite 
wave problems in the case of large wave number (cf. [2 2) and the references therein) . 

This paper is an attempt to address the first difficulty mentioned above for the 
Maxwell equations. In particular, our goal is to design and analyze discretization 
methods which have superior stability properties and give optimal rates of convergence 
for the Maxwell problem. Motivated by our previous experiences with the Helmholtz 
equation [10[ 111) , we again try to accomplish the goal by developing some interior 
penalty discontinuous Galerkin method for problem (|l.l|) - ()1.2p . The focus of the 
paper is to establish the rigorous stability and error analysis for the proposed IPDG 
method, in particular, in the preasymptotic regime (i.e., when k^h > 1). For the ease 
of presentation and to better present ideas, we confine ourselves to only consider the 
linear element in this paper and will discuss its high order extensions in a forthcoming 
paper. 

The remainder of this paper is organized as follows, section [5] is devoted to 
the study of the coercivity of the Maxwell operator and the wave-number explicit 
estimates for the solution of (|1.1[) - (|1.2[) . We show that the sesquilinear form asso- 
ciated with the Maxwell problem satisfies a generalized weak coercivity (i.e., inf-sup 
condition). This coercivity in turn readily infers the wave- number explicit solution 
estimates which were proved in [8l[T3]. We note that the proofs of both results given in 
this paper are of independent interest and refer the reader to 9 for further discussions 
in the direction, section |3] presents the construction of our IPDG method and some 
simple properties of the proposed discrete sesquilinear form, section 3] studies the co- 
ercivity of the discrete sesquilinear form and derives stability estimates for the IPDG 
solutions. It is proved that the discrete sesquilinear form satisfies a coercivity for all 
mesh size h > and all wave number A; > and for general domains including non- 
star-shaped ones, which is stronger than the generalized weak coercivity satisfied by 
its continuous counterpart. All these are possible because of the special design of the 
discrete sesquilinear form and the special property curl curl v/i = (element- wise) 
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for all piecewise linear functions v^. This coercivity in turn readily infers the well- 
posedness and stability estimates for the discrete problem without imposing any mesh 
constraint, section [5] devotes to the error analysis for the proposed IPDG method. 
By using the discrete stability estimates and adapting a nonstandard error estimate 
technique of [lOj . we derive both the energy-norm and the L^-norm error estimates 
for the IPDG method in all mesh parameter regimes including pre- asymptotic regime 
(i.e., k^h > 1). Finally, we present some numerical experiment results in section [6] 
to gauge the theoretical results and to numerically examine the pollution effect (with 
respect to k) in the error bounds. 

2. Generalized inf-sup condition and stability estimates for PDE so- 
lutions. The standard space, norm and inner product notation arc adopted in this 
paper. Their definitions can be found in [2111]. In particular, (•, •)q and (•, •)s for 
Q n and S c dfl denote the L^-inner product on complex-valued L^{Q) and i^(S) 
spaces, respectively. For a given function space W, let W = (W)^ . In particular, 
L2(rj) = (L2(f^))3 and H'=(rj) = {H''{n)f. We also define 

H(curl, n) := {v e L^{n); curl v e L^in)}, 
H(div,r2) := {v6 L2(r2); divve L^{n)}, 
U{divo,n) ■■= {v6 L2(1]); divv = 0}, 

V := {v e H(curl, 12); vt e L^{T)}, 

V := {v e H(curl, fl) ; curl v e H(curl, n),ve H(curl, F) } . 

Throughout this paper, the bold face letters are used to denote three-dimensional 
vectors or vector- valued functions, and C is used to denote a generic positive constant 
which is independent of h and k. We also use the shorthand notation A < B and 
B > A for the inequality A ^ CB and B 5= CA. A ~ B is a. shorthand notation for 
the statement A < B and B < A. 

We now recall the definition of star-shaped domains. 

Definition 2.1. Q cz R'^ is said to be a star-shaped domain with respect to 
€ Q if there exists a nonnegative constant cq such that 

(2.1) {ji-jiQ)-UQ^CQ Vxe(5g. 

Q c R"^ is said to be strictly star-shaped if cq is positive. Where vq denotes the 
unit outward normal to dQ. Throughout this paper, we assume that is a strictly 
star-shaped domain. 

Introduce the following sesquilinear form on V x V 

(2.2) a(u, v) := (curl u, curl v)^ - fc^(u, v)j2 - iA<UT,VT>r, 

Then the weak formulation for the Maxwell system (|l.ip - (|1.2p is defined as seeking 
E e V such that 

(2.3) a(E,v) = (f,v)j^ +<g,VT>r Vv e V. 

Using the Fredholm Alternative Principle it can be shown that problem (|2.3|) has a 
unique solution (cf. 0|T7]). 

Note that choosing v = VV' with ip e -ffg (51) shows that (fc^E -I- f , V'0)n = 0, or 



(2.4) 



div(fc2E + f ) = inn. 
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Next, we prove that the sesquilinear form a(-, •) satisfies a generalized weak coer- 
civity which is expressed in terms of a generalized inf-sup condition. 

Theorem 2.2. Let n cz be a bounded star-shaped domain with the positive 
constant cq and the diameter R = dim(r2). Then for any u e V n H(divo,il) there 
holds the following generalized inf-sup condition for the sesquilinear form a(-, •).• 

. ^, llmafu, v)| |Rea(u, v)| 1., 

(2.5) sup J II I + sup ,,, 5? -Me, 

where 

(2.6) 7 := max{4/fci?, M}, M := \ \ 

Acq 

(2.7) |||u|||L2(a) := (fc'HuH^ .^^^ + ecnMl.^^^^y, 

1 

(2.8) \\u\\e := (fc'||u||^2(jj) + k^cnMl^f^T) + il curlujl^.^^^ + co|| curl ujl^.^j,^) ' . 

Proof. Let w := x — xji- Setting v = u in (|2.2p and taking the real and imaginary 
parts we get 



(2.9) Rea(u,u) = || curlujl^.^^^ - fc2[|^j|2^^^^^ 

llL2(r)- 



(2.10) Ima(u,u) = -A||uTp 



Alternatively, setting v = curl u x w in (|2.2p (notice that v e V is a valid test function 
for u e V), taking the real part, and using the following integral identity (cf. [5]) 

(2.11) Re(u, v)n + il|ul|2,(^,j + l<w • IV, |u|2>r 

= Re(w X u, u X i/)r + Rc(div u, u • w)o 
and the assumption that div u = 0, we get 

(2.12) 2 Re a(u, v) = 2 Re( curl u, curl v)o - 2k'^ Re(u, v)n + 2A Im<UT, VT>r 

= 2Re(curlu, curlv)o + A:2||u||l2(q) + fc2<w • iv, |up>r 
— 2k'^ Re(w X u, u X iy)r + 2A Im(uT, VT)r- 
From (|2.9p and (|2.12p and using the following integral identity (cf. [5]) 

(2.13) 2 Re( curl u, curl v)o = 1| curl u\\12(^q) + <w • | curl up>r, 
we have 

(2.14) 2/c2||u||2,(^) = k^uWl^^^ + k^uWl^^) 

= -2Re(curlu, curlv)n - k'^(-w ■ v, |up>r 
+ 2fc2 Re(w X u, u X u\ — 2 A Im(uT, VT>r 
+ 2Rea(u, v) + || curlu||L2(Q) — Rea(u, u) 

= -(w • i>, I curlup>r - fc^<w • v, |up>r 

+ 2fc2 Re(w X u, u X u)t — 2 A Im(uT, VT)r 
+ 2 Re a(u, v) — Re a(u, u) 
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= -<w • u, I curl up>r - /i;^<w • v>, |up>r 

- 2fc^<w ■ v,\ux i^P>r + 2A:^ Re<WT x u, u x v)r 

— 2AIm(uT, VT)r + 2Rea(u, v) — Rea(u, u). 

Here we have used the decomposition w = (w • + to obtain the last equahty. 

On noting that (wt x u, u x v>}-p = (u • i/,wt ■ UT}r, \\'^\\L'^(n) ^ and that 
|v| ^ |curlu||w|, using the star-shaped domain assumption and Schwarz inequahty 
we obtain 

(2.15) 2k^\\u\\l2^nj =S -coll curlu||2,(j,j - /c2^n||u||^2(r) 

-2fc2cn||uT||^2(r) + 2fc2i?|| 

u||L2(r)l|uT||L2(r) 
+ 2Ai?||uT||L2(r)ll curlu||L2(r) + 2Rea(u,v) — Rea(u,u). 

curlu||^,(,) - ^||u||^,(,^ - 2A;2co||uT||L(r) 

H l|uT||L2|'r'i + 2Rea(u,v) — Rea(u,u). 

cn ^ ' 

Finally, it follows from ([Q}]) and (|2l^ that 

(2.16) 2fc2||u||L2(n) + 2|| curlu||L2(n) + co|| curlu||^2(r) + fc^cn||u||^2(r) 
^ M|Ima(u, u)| + |Rea(u,4v)|, 

where v = curl u x w and M is defined in (|2.6p . 

It is easy to check that there holds for v = curl u x w 

|||v|||L2(fj) fcii-llulU. 

Hence, it follows from (|2.16p that 

/2 |Ima(u, u)| _^ |Rea(u,4v)| ^ |ImQ(u,u)| _^ \ Rea(u,4v)| 
Me II|4v|||l2(o) " ||u||_E 4/ci?||u||£; 

1 Af|Ima(u,u)| + |Rea(u,4v)| 1 

7 llullij 7 

where 7 = max{4fci?,M} as defined in (|2.6p . The proof is complete. □ 

An immediate consequence of the above generalized inf-sup condition is the fol- 
lowing stability estimate for solutions of problem (|l.ip - (|1.2l) . 

Theorem 2.3. In addition to the assumptions of Theorem \2.S[ assume that 
f e H(div, n) and g e L^(r). Let E G V n H(div, f2) be a solution of the variational 
problem (j2.3|) . Then there holds following stability estimate: 



(2.18) II curlE||L2(f2) + A:||E||L2(f2) + Vcoll curlE||L2(r) + fcv^l|E||L2(r) 

<fc-V^(f,g) + fc"'ll divf||L2(o) 

for all k,\ > 0. Where 

(2.19) M(f,g) := ||f||L2(f2) + c-^||g||L2(r). 
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Proof. Let ip e Hq (il) solve 

(2.20) A(p = -k-^divi inn. 

Set F = V(/3 and u = E — F, where E is a solution to (|2.3p . Trivially, we have 
curl F = and div F = -k~^ div f in fl, and Ft = Vt"^ = on T. By we also 
have divu = div(E — F) = 0. Hence, u e H(divo,ri). Moreover, since E satisfies 
(|2.3p . it is easy to verify that u satisfies 

(2.21) a(u,v) = (f + fc2F,v)o + <g,VT>r Vv e V. 

Testing (|2.20[) by tf and integrating by parts on both sides of the resulting equation 
yield 

WvWhm = -fc"'(f, V(^)j2 ^ A;-2[|f[|L2(o)[|V(/p[|L2(a). 

Hence, 

(2.22) l|F!|L^(n) = ||V(^||L2(n) < k-^r\\j^2^n)- 

Alternatively, testing (I2.20p by Vip • w = F • w with w = x — xjj , using the following 
Rellich identity for the Laplacian (cf. pUl \7\): 

2Re(A(pV^-w) = |V<^p + 2Re(div(V(^ • w)) -div(w|V<^n, 

and integrating by parts we get (note that F^ = 0) 

-2fc"2(div f , F • w)o = Ml^t^a) + 2 Re<F • iv, F • w>r - <w • \F\'^}r 
= ||F||^,(^) + <w«.,|F|2>r. 

Hence, by (I2.22p and the star-shaped domain assumption we obtain 

(2.23) mhin) + cn\\F\\l.^r) ^ 2fc-2||w|Ucc.(n)|| divf ||l.(o)||F||l2(o) 

sS 2k^^R\\ divf||L2(f2)l|fl|L2(f2)- 

Finally, by p.22p and Schwarz inequality we get 

(2.24) |(f + fc^F, v)o + <g, VT>r| < 2fc-iM(f , g) [k^w\\l,^^^ + fc'co||vT||L(r)) ' 

^2fc-iM(f,g) |||v|l|L2(n). 
It follows from the generalized inf-sup condition (12.51) . (|2.2ip and (|2.24p that 

(2.25) 7-i||u||£ «;4fc-iM(f,g), 

which together with ()2.23p and the relation u = E — F as well as the definition of the 
energy norm Hujl^; infer that (again, note that Ft = 0) 

||E|U ^ ||u|U + ||F||b 

4fc-i7Af(f,g) + fc(||F||2,(^^ + co||F||2,(j,))^ 
< 4fc-i7M(f,g) + 2i?t|f||L2(o) + (2A;)-2 [| divf||L2(n). 
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Hence, (I2.18|) holds. The proof is complete. □ 
We conclude this section with a few remarks. 

Remark 2.1. Since problem (11.11) - (jl.2[) is linear, the stability estimate (j2.18l) 
immediately implies the uniqueness of the problem in the function class in which the 
estimate is derived. This provides an alternative method (to the traditional integral 
equation method and the unique continuation method) for establishing uniqueness (and 
existence) for the Maxwell problem (II. 1[) - (11.21) . 

Remark 2.2. (a) The generalized inf-sup condition (|2.5p is a stronger result than 
a stability estimate for the solution of the Maxwell problem. The reason to restrict 
u e H(divo, 51) in (|2.5p is that curl operator has a non-trivial kernel. 

(b) Stability estimates similar to (I2.18[) were established independently early in 
0/ and U3^ . (|2.18|) also explicitly shows the dependence on the size and the shape 
constant of the domain. Such an estimate plays an important role for designing 
multilevel Schwarz preconditioners for discretizations of (|2.3p and for doing practical 
simulations because in practice the size of the computational domain Q is often taken 
to be proportional to the wave length. 

In addition, not only the sharp wave number- explicit and domain size- explicit 
stability estimate (j2.18p is obtained as a corollary of the generalized inf-sup condition 
(j2.5p . but also the derivation reveals some deep insights about the dependence of the 
solution on the datum functions and the domain. 

(c) The generalized inf-sup condition (|2.5p provides a guideline for constructing 
"good" numerical schemes for the Maxwell equations. We shall call a discretization 
method "a coercivity preserving method" if it satisfies a discrete inf-sup condition 
which mimics the continuous inf-sup condition. Constructing such a coercivity pre- 
serving IPDG method is one of primary goals of this paper. 

(d) Generalized inf-sup conditions similar to (j2.5l) also hold for the scalar Helmholtz 
equation and the elastic Helmholtz equations (cf. ^). 

Based on the above stability estimates in lower norms, one can also derive stability 
estimates in higher norms when the solution E is sufficient regular. We state an H^- 
estimate for curlE below without giving a proof (cf. [13l Remark 4.9]). 

Theorem 2.4. Suppose that divf = and the solution E of problem (|l.ip - (|1.2p 
satisfies E e H*(curl, il) for ^ < 6 ^ 1. Then there holds estimate 

(2.26) l|E||Hncuri,o) ^ (l + A + fc)Af(f,g) + ||g||^i^^^, 
where 

(2.27) H''(curl, ft) := {u e H*(rJ); curl u e H'^(17)}, 

1 

(2-28) l|u||H'5(curi,n) := (hllH-'in) + II curlu||^,(s^)^ ' . 

3. Formulation of discontinuous Galerkin methods. To formulate our IPDG 
methods, we first need to introduce some notation. Let {Th} be a family of partitions 
(into tetrahedrons and/or parallelepipeds) of the domain fl parameterized hy h > 0. 
For any "element" K e Th, we define hx ■= dia.m{K). Similarly, for each face T 
oi K e Th, define hjr := diam(J^). We assume that the elements of Th satisfy the 
minimal angle condition. Let 

Sl := set of all interior faces of Th, 

:= set of all boundary faces of 7^ on F = d^. 
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We define the jump [v] and average {v} of v on an interior face T = dK n dK' as 

r-,-|| [ v|x - v|i^-, if the global label of K is bigger, , , , 1 , , , x 

1 v|a" -v|k, if the global label of if' is bigger, ^^'"^ ' 2*^'^ 

If J" e 5^, set [v]|jr = v|jr and {v}|^ = v|jr. For every T = (5if n dK' e let t'j^ 
be the unit outward normal to the face T of the element K if the global label of K is 
bigger and of the element K' if the other way around. For every T e , let i/jf = v> 
the unit outward normal to dVl. 

To formulate our IPDG methods, we recall the following (local) integration by 
parts formula: 

(3.1) (curl E, F)k = (E, curl F) ^ - <E x i^^, Ft}8k- 

where Ft = {i'k x F) x v>k- 

Next, multiplying equation (jl.ip by a test function F, integrating over K e Th, 
using the integration by parts formula p.ip , and summing the resulted equation over 
all K e Th we get 

(3.2) 2 ((curlE,curlF)K-<curlExiyA,FTW)-fc2(E,F)n = (f,F)a. 

To deal with the boundary terms in the big sum, we appeal to the following algebraic 
identity. For each interior face J- = K n K' e there holds 

(3.3) <curlE X iv^, FtV + <curl E x uk',Ft}j^ 

= <[curlE X u^l {Ft}>^ + <{curlE x j^^}, [Ft]>^. 

Substituting identity p.3p into (j3.2p after dropping the first term on the right-hand 
side of p.3p (because [curlE x I'jrJIjr = if E is sufficiently regular) yields 

2 (curl E, curl F) A' - ^ <{curl E x z^^}, [Ft]>^ 

- < curl E X I., FT>r - ^^(E, F)n = (f , F)o, 

Utilizing the boundary condition (|1.2p in the third term on the left-hand side and 
adding a "symmetrization" term then lead to the following equation: 

(3.4) 2 (curl E, curl F) A - ^ (<{curl E x z^^}, [Ft]>^ 

Ken resi 

+e<[ET],{curlF x «^^}>^) - iA<ET, FtV - fc'(E, F)^ = (f , F)f^ + <g, FtV 
where e = —1, 0, 1. 

The most important and tricky issue for designing an IPDG method is how to 
introduce suitable interior penalty term(s) on the left-hand side of (13. 4p . Obviously, 
different interior penalty terms will result in different numerical methods. As it was 
proved in fT5^, using the standard interior penalty terms will lead to IPDG methods 
which require a restrictive mesh constraint to ensure the stability and accuracy in the 
case of large wave number k. Inspired by our previous work |10| on IPDG methods for 
the Helmholtz equation and guided by our stability analysis (see sectional), here we 
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introduce some non-standard interior penalty terms into (j3.4l) , which we shah describe 
below, and the IPDG method so constructed will be proved to be absolutely stable 
(with respect to wave number k and mesh size h) in the next section. 

To define our IPDG methods, we first introduce the "energy" space V and the 
sesquilinear form 6^(-, •) on V x V as follows. 

V := n ■■= {v e H(curl,/^); ^\sk e I^H^K), curW\sK 6 L^{dK)}. 

KETh 

(3.5) bl{u,v) := 2 (curl u, curl v) if 



(3.6) Jo(u,v):= 2 ^<[ut],[vt] 



(^({curlu X vjr}, [vt])^ + e([uT], {curlv x i^jr})^ 
i(Jo(u,v) + J'i(u,v)), 

(3.7) Ji(u, v) := 2 7i,^ft.;r([curlu x n^], [curlv x n^]>^, 

where 7o,j^ and 71. j^- are nonnegative numbers to be specified later. 

Remark 3.1. (a) Clearly, 6^(-, •) is a consistent discretization for curl curl since 
(curl curl u, v)n = b\{u, v) for all u e TrP{il) and v e V with v^lr = 0. 

(b) The terms in — i(^o(u, v) + J7i(u, v)) are called penalty terms. The penalty 
parameters — i7o,j^ and — i7i,j^ are pure imaginary numbers with negative imaginary 
parts. Our analysis still applies if they are taken as complex numbers of negative 
imaginary parts. 

(c) The Jo term penalizes the jumps of the vector field u and the J\ term penalizes 
the jumps of the tangential component of the vector field curlu. which, to the best 
of our knowledge, has not been used before in the context of IPDG methods for the 
Maxwell equations. They play a vital role for our IPDG methods being absolutely 
stable, see section^ 

(d) e = —1,0,1 correspond to the nonsymmetric, incomplete, and symmetric 
IPDG methods for the Poisson problem. In the remainder of this paper, we shall 
only consider the symmetric case e = 1 and set bh{-, •) = b\{-, ■) for notation brevity. 

With the help of the sesquilinear form 6/1 (•, •) we now introduce the following weak 
formulation for (frT|) - (fr2|) : Find E e V n H(curl, Q) such that 

(3.8) a,,.(E,F) = (/,F)o + <g,FT>r VF e V n H(curl, 1]), 
where 

(3.9) a,,(E, F) := &,,(E, F) - fc2(E, F)o - iA<ET, FtV- 

From (123), it is clear that, if E e is the solution of (ITTT]) p^ . then (ISTSI) 

holds for all F e V. 

For any K e 7ft,, let Pr{K) denote the set of all complex- valued polynomials whose 
degrees in all variables (total degrees) do not exceed r(> 1). We define our IPDG 
approximation space Vft as 

Vft := Yl Pr{K). 



10 



X. FENG AND H. WU 



Clearly, V,^ c V c L^{n). But V,,. cj: H(curl,17). 

We are now ready to define our IPDG methods based on the weak formulation 
((3:8)) : Find Eh e Vh such that for all Fh e V,^ 



(3.10) ah(Eh,Fh) = (/, Fh)n + <g, (Fh)T\. 

We note that p.lOp defines a family of IPDG methods for r 5= 1. For the ease 
of presentation and to better present ideas, in the rest of this paper we only consider 
the case r = 1, the linear element case. In the next two sections, we shall study the 
stability and error estimates for the above IPDG method with r = 1. Especially, 
we are interested in knowing how the stability constants and error constants depend 
on the wave number k (and mesh size h, of course) and what are the "optimal" 
relationship between mesh size h and the wave number k. We remark that the IPDG 
method with r = 1 uses piecewise linear polynomials even for Cartesian meshes. By 
contrast, for the corresponding linear conforming edge element method on Cartesian 
meshes, the trial functions have to be chosen as piecewise trilinear polynomials. We 
also note that the linear system resulted from (I3.10p is ill-conditioned and strongly 
indefinite because the coefficient matrix has many eigenvalues with very large negative 
real parts. Solving such a large linear system is another challenging problem associated 
with time harmonic Maxwell problems, which will be addressed in a future work. 

For further analysis we introduce the following semi- norms/norms on V: 

(3.11) llcurlvll^^^.^^^ := ^ || curl vl^^^^^, 

KeTh 

(3.12) : = ||curlv||2,(^^^ + ||v||2,(^^ 

'^o.-^ii r„ 1 ii2 , „ L 11 r„„„i„ ., 1 ii2 



+ Zl (x^ll [^^] + ^i.-^^-^ll [curlv X u^] 

= 11 curlv||2,(.^^) + ||v||2,(^) + Jo(v,v) + Ji(v,v), 

(3.13) M\l^ := j|vf^^ + ^ {curlv x «.^} ||2 

Clearly, the sesquilinear form bfi{-, •) satisfies: For any v e V 

(3.14) Re6„(v,v) = 1| curl vjl^.^^^j - 2 Re ^ <{curlv x z^^},[vt]>^, 

(3.15) Im6;,(v,v) = -Jo(v,v) - J'i(v,v). 



4. Discrete coercivity and stability estimates. In this section we shall prove 
that the discrete sesquilinear form a/i(-,-) satisfies a discrete coercivity, which is 
slightly stronger than the generalized inf-sup condition proved in the previous sec- 
tion for the sesquilinear form a(-,-). Such a discrete coercivity is possible for the 
linear element because curl curl v^j = (defined element- wise) for all G Vh- As 
an immediate corollary of the discrete coercivity, we shall derive a priori estimates for 
solutions of p.lOp for all h,k > Q, which then infer the well-posedness of (I3.10p . 

We state the first main theorem of this section which establishes a coercivity for 
the discrete sesquilinear form ah{-, •)■ 

Theorem 4.1. Let 70 = min^r^^i {70,7^}, 7^ = min^r^^/ {71,;^}, and /imin = 
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minjTgg^ {hjr}. Then there exists a constant < C < 1 such that 

C 

(4.1) \ah{uh,Uh)\^ —WuhWl h Vu^eV,, 

for all A:,7o,7i > 0. Where 

. , 1 11 

:= YT, + 1.2U2 + — + 1. 

A/lmin 7l« /i^in 70 

(4.3) \\uh\\E,h := (ll curluft||^2(r,) + A:'||uft||L2(n) 

+ Jh{jQ{Uh,Uh) + Ji{Uh,Uh) + Al|(u;i)T||L2(r) 

Proof. For any J" e S^, define := (J e 7^; cKnF ^ 0}. By (EH), ([XTi)) . 
and the following trace inequality 

(4.4) ||K}||l2(^) «:CV^||v„||L2(n^) Vv„ e 
for some /ijr-independent positive constant C, we get 

(4.5) Rea,,(u,„u,i) || curl u?,||^2(7-^) - fc2||u,,||^2(f2) 

+ 2 2 ||{curlu/, X /y^}I|L2(^)||[(Uft)T]||L2(^) 
+ ^2 V' II CUrlUft||L2(o^)||[(u,0T]||L2(^) 

=5 -||curlu,,||2,^^^j -A:2||^^||2^^^^ 2 II [(u'«)t] IIl2(^)- 

Since U;j e V/j is piecewise linear, then curl curl U/j = in each K e Th- By 
integrating by parts and using the trace inequality (14.41) we obtain 



curluft||L2(7^^) = 2 (curl u/i, curl Uft)^ 



Yj ((curl curl u/i,u,i)^ - (curl u,, x VK,{^h)T)g 



KeTh 



dK 



= - 2 (curlu,, X VK,{nh)T)j, 

- Yj (<[curlu,j X iy^],{(u,,)T}>_^ + ({curlu,, X iy^},[(u/,)T]>^ 

Yj V' II <^"''l"'»llL2(n^)ll("'»)7^llL2(.7^) 
+ '^1] II[CUI"1U,, X Zy;r]l|L2(;r)l|u,,||L2(n^) 

+ Zl ^^'l|curlU;i||L2(o^)||[(u,,)T]||L2(^) 
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^-\\curluh\\l2^r^) + —\\uh\\l2(^<^^ + C 2] V^IK)tIIL(^) 



Hence, 

^2 

(4.6) 2||curlu^||2,(^^) <; — ||u^||2,^^^ +C ^ /i^'IIK)TllL(^) 



Adding (I4.5P and (I4.6P and rearranging the terms yield 



(4.7) ||curlu,,||2,(^^) + fc2||u,^j|2^ 



in) 



=g-2Rea^(u^,u,.) + -^A||(u^)T||L(r) + - ^ ^11 [(u^M 



c 



+ ^ p/,2 E 71,^ VII [curium X i/^] ||2 

Therefore, by the definitions of J'o(-, •) and ^i(-, •) and the identity p. 151) we get 

II curlufi||L2(-7^j + fc2||uft||L2(n) + 7ft(Jo(uh, u,i) + Ji(uft,u,i) + A||(u,i)T||L2(r)) 
^ -2Reah{uh,Uh) + Cjh{Jo{uh,Uh) + Ji{uh,Uh) + A||(u,i)T||L2(r)) 
= -2Reah{uh,Uh) - C7/1 Im a/i(u/i, u^) 

< C7/i(|Rea,i(u,i,u/i)| + | Iniaft(u,,, u/i)|) Cjh\ah{uh,Uh)\, 

where 7^ is defined by (|4.2p . Hence, (j4.ip holds. The proof is completed. □ 

Remark 4.1. ('aj The discrete sesquilinear form ah{-,-) satisfies a stronger co- 
ercivity than its continuous counterpart a{- , ■) does, see Theorem \2.2i Moreover, the 
proof of Theorem \4-l\ is simpler than that of Theorem \2.2\ all these are possible be- 
cause of the special form of a/i(-, •) and the fact that curl curl v/j = in K e Th for 
all piecewise linear functions V/j G . However, a weak coercivity is only expected to 
hold in the case of high order elements. 

(b) It is also important to point out that Theorem \4. 1\ holds without assuming that 
is a star-shaped domain. 

An immediate consequence of the above discrete coercivity are the following a 
priori estimates for solutions to the IPDG method p.lOp . 

Theorem 4.2. Every solution 'Eh of the IPDG method p.lOp satisfies the fol- 
lowing stability estimates. 

II curlE/i||L2(r^) + fc||Eft||L2(j2) < fc~^7/i||f||L2(o) + (A^S/0^ l|g||L2(r), 
{Jo{-Eh,-Eh) + Ji{Eh,Eh) + A||(Eft)T||^2(r))^ < fc^Sj l|f IIl2(o) + A~^||g||L2(r). 
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Proof. By (j3.10[) and Schwarz inequality we get 

|a,,(E,„E,i)| = |(f,E,,)f2 + <g, (Eft)T>r| *S \\thHn)\\Ehh^n) + l|g||L2(r)||(E/i)T||L2(r) 

(fc'l|E,.ilL(o) + A7^i|(E,0|||L^(r))Mfc"1f|lL^(o) + (A7h)-'||g|lL(r))^ 
llE^llij,,, (fc"^||f||L2(n) + (A7/i)"^l|g[|L2(r))- 



The desired estimates follow from combining the above inequality with (|4.ip . The 
proof is completed. □ 

The above discrete stability estimates in turn immediately imply the well-posedness 
of the IPDG method (|XTU)) . 

Corollary 4.3. There exists a unique solution to (|3.10p for any fixed set of 
parameters k, hjr^ 7o,J^7 li.J" > 0. 

5. Error estimates. In what follows, we suppose 70, jr — 70 and 7i_jr ^ for 
brevity. For simplicity, we assume that div f = and that Th is a quasi-uniform 
partition of consisting of tetrahedrons. Let h := max{hj^; K e Th}- 

5.1. H(curl, r2)-elliptic projection and its error estimates. Let E be the 

solution to problem (|l.ip - (|1.2p and Eh e V/i be its IPDG H(curl, r2)-clliptic projec- 
tion defined as follows. 

(5.1) bh{E-Eh,Vh) + {'E-Eh,Vh)n = Vv,, e V,,. 

The following lemma establishes the continuity and coercivity for the discrete 
sesquilinear form bh{-,-). 

Lemma 5.1. For any v,w e V, the sesquilinear form hh{-, ■) satisfies 

(5.2) |6,i(v, vif) + (v, w)o|, |6/i(w, v) + (vif, v)j2| < |||v|||ug|||w|1|£)g. 
In addition, there exists a positive constant 7 such that, for 70 5= 7, 

(5.3) Rebh{-Vh,Vh) - Imbh{vh,vh) + {vh,vh)n ^ ^\\\vh\\\DG Vv/i e V^. 

Proof Clearly, ([53) follows from the definitions ([331) -([SJl), ((3lT1) - ((XT3| . and 
Schwarz inequality. It remains to prove (|5.3p . 
From (|5T1|) - (I5T^ we have 

Rebhivh,Vh) -lmbh{vh,Vh) + (v/i,v/,.)o 

= |||v;,|l|2,^-2Re <{curlv;,xi.^},[(v;,)T]>, 

- y. —\\ {curlvft X ujr} [|2 

It follows from the derivation of (|4.5p that there exists a constant cq > such that 
2Re 2 <{curlv^ x u^} ,[{vh)T]\ 
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On the other hand, from (|4.4p . there exists a constant ci > such that, 

^6£^ 70,e 70 

Therefore, 

Rebh{vh,Vh) - lmbh{-Vh,Vh) + (v,i, v/Jo 5= f 1 - 7 - ^" '^M |||vh||||,G 

V 4 70 / 

which gives (|5.3p if 70 large enough. The proof is completed. □ 

Remark 5.1. The coercivity and continuity of bh{-,-) ensure that the above 
Uncurl, D,)- elliptic projection is well defined. 

The following lemma establishes error estimates for E — Eh- 
Lemma 5.2. Suppose problem (jl.ll) - (|1.2p is H"^ -regular, then, under the condi- 
tions of Lemma \5. 11 there hold the following estimates: 

(5.4) |||E-E,,|||i3G < /^(l + 7l)^l|E||Hi(curl,0), 

(5.5) ||E - E^||L2(n) < h^il + 7i)7e(E), 

(5.6) ||E - E,,||L2(r) <h^l + ji)TZ{E), 

where 

(5.7) 7^(E) := (l + 7i)^l|E||Hi(curi,o) + MnHn)- 

Proof Step 1: It follows from [HI [HI HI] that there exists E,i e V,, n H(curl, Vt) 
(i.e., the conforming Nedelec interpolation of E) such that the following estimates 
hold: 

(5.8) ||E-E,,||L2(n) < /i'||E||ff2(n), 

(5.9) ||E-E,,||L2(r) < /i^||E||ff2(o), 

(5.10) ||E - E;i||H(curl,0) /l||E||Hi(curl,a) j 

(5.11) |||E-E„|||cG </»(!+ 71)^ l|E||Hi(curl,0), 

where l|5.1ip can be proved by (j5.10p . the commuting property between the curl- 
conforming interpolation operator and the div-conforming interpolation operator [171 
Lemma 8.13], and the trace inequality. 

Let := E/i - E,, and := E - E,,, then E - E,, = 'S'/i - By (|54|) we 
have 

(5.12) bh{^h, */0 + (*/., */.)a = 6/.(*h, ^h) + (*h, ^h)n. 
Step 2: From Lemma 15.11 and (|5.12l) we get 

(5.13) i|||*^|||2,c <Re6,(*,,*,0 - Im 6, (*,,*„) + ||*,||^,(^) 

= Re(6,i(*,i, *,,) + (*,,., *,0o) - lT[v{bh{^h,^h) + (*/i, ^h)n) 
<|||*,,.|||dg |||*/i|||_DG. 
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Therefore, it foUows from (|5.1ip that 

(5.14) lll*h|l|DG < l||*/>.|bG < /i(l + 7l)'l|E||Hi(curl,0), 

which together with the relation E — Eh = "i^h and the triangle inequality 

immediately infer (j5.4p . 

Step 3: To show (|5.5p . we first need the following results that can be proved 
by following the proof of jl5| Proposition 4.5] and their proofs are omitted: for any 
Vh e Yh there exists e n H(curl, ft) such that 

(5.15) ||v„ - v^J|^,(,^) Y VII [(v„)t] IIL(^), 

(5.16) II curl(v, - ^DWh^r,) Y h^'W [K)t] Wl^^y 

Let e V;i n H(curl, ft) be the conforming approximation of as defined above. 
Then it follows from the definition of the norm \\-\\]jq (cf. (I3.12p ). the above two 
estimates, and (I5.14p that 

(5.17) II*, - mi^Hu) + h\\ curl(*, - *DllL^(r.) ^ %^h\\^h\\DG ^ h'^^). 
Noting that 

||E - EhWl.^^^ = (E - E,„ E - E,.)^^ - (E - E,, *^)^^ - (E - E,, - *^,)^, 
we have 



(E-E,^)^ 

l|E — Eh 



l|E - E,i||L2(f2) ^ [|E - Eft,||L2(j2) - — ^— h l|*/i - *^||L2(f2), 



which together with (|5.8p and (|5.17p yields 
(5.18) ||E-Eft||L2(n) </i'7^(E) 



2^,^, lE-E,,*^)^ ^ 

II E - E/i||L2(f2) 



Step 4-' We need to bound the last term on the right-hand side of (|5.18p . Notice 
that e V/i n H(curl, $7), by using a standard duality argument, see Appendix, 
based on the Helmholtz decomposition of we can show that 

(5.19) ~ \ " "'^ < (1 + 7i)/i'7l(E). 

I|E - E/i||L2(f2) 

Step 5: The desired estimate (|5.5p follows from combing (|5.18p and (|5.19p . Finally, 
dill) follows from ||E - E?jL2(r) ||E - Eh||L2(r) + ||E,, - E/jL2(r), (EH), the trace 
inequality, (|5.8p . and (|5.5p . The proof is complete. □ 

Remark 5.2. T/ie Pi -conforming Nedelec edge element (of second type) projec- 
tion E/i of E is introduced and used in the proof to simplify the analysis at the expense 
of requiring Th to he a quasi-uniform and conforming mesh. We note that the proof 
is still valid if one replaces the Pi-conforming Nedelec edge element projection by the 
Pi-IPDG projection without assuming Th is a quasi-uniform or conforming mesh. As 
expected, the new proof will he more complicated and technical, and is left for the 
interested reader to explore. 
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5.2. Error estimates for IPDG method (|3.10l) . The goal of this subsection 
is to derive error estimates for scheme p.lOp . Instead of using the well-known Schatz 
argument [H] [TSl [Ml [H] , which is the (only) technique of choice for deriving error 
estimates for indefinite problems in the literature, we shall obtain our error estimates 
by exploiting the linearity of the Maxwell equations and making strong use of the 
discrete stability estimates proved in Theorem 14.21 and the projection error estimates 
established in Lemma 15.21 This new technique, which is adapted from [lOj . allows us 
to derive error estimates for e/j := E — Eh without imposing any mesh constraint. 

It is easy to check that there holds the following error equation: 

aft(e/i, v,i) = Vv/jeV/i. 
Let r/^j := E — Eh and ^h ■= E/i — E/j, then eh = Vh^ ^h- From (|5.ip we get 

(5.20) ahi^h^'^h) = ah{'nh,^h) = bhiVh^'^h) - k^{r]h, ^^h)n - iX{iVh)T, (vh)T>r 

= -(fc2 + l){r]h, vh)n - iA<(?7^)T, (v,,)T>r Vv^ e V^, 

The above equation implies that ^h ^ is the solution of scheme p.lOp with the 
source functions f = — (fc^ + l)rih and g = ~X{rih)T- Hence, an application of 
Theorem 14.21 and Lemma [5^ immediately infers the following estimate for 
Lemma 5.3. Under the conditions of Lemma \5. 1[ there holds 

(5.21) \\UDG + mhhHn)^CsUl + ^i){k^h^ + Xhl)niE), 
where 

(5.22) Csta := max (fc-i(l + 7,,), {X-\l + lh))~'), 

and 7/1 is defined by (|4.2p . 

By Lemmas 15.21 and 15.31 and the triangle inequality we then obtain the following 
main theorem of this section. 

Theorem 5.4. Let E and E/i be the solutions to problem (ll.ip - ()1.2|) and scheme 
()3.10|) . respectively. Assume E e H^(ri). Then, under the conditions of Lemma \5.1\ 
there hold the following error estimates: 

(5.23) ||E-E,,||i5G < {h + C,Ul + li){k'h^ + Xhi))n{E), 

(5.24) ||E-E,,||l2(o) < (/^2 + 5,tafc-l(fc2/^2 + A/^i))(l + 7l)7^(E). 

To bound 7?.(E) in terms of the source functions f and g, we need to bound 
l|E||H2(n) and l|E||Hi(curi,n) by the source functions. To the end, we appeal to the 
solution estimate (|2.26p to get 

(5.25) n{E) < (A + fc)Af(f,g) + l|g||^i(^) 

Substituting (|5.25p into (I5.23P and (|5.24p yields the following explicit in all pa- 
rameter error bounds for E — E;j. 

Corollary 5.5. Suppose fc, A > 1, and < 71 < 1. Under the assumptions of 
Theorem \5.4\ there exist constants Ci and C2 independent of k,X, and h such that 

(5.26) ||E - EhWoG ^ Ci{k + X)h + C2Csta(fc + X){k^h^ + Xhi), 

(5.27) ||E-E^||l2(j2) < Ci{k + X)h^ + C2Cst.k-\k + X){k^h'' + Xhi). 
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Remark 5.3. (a) If X = 0{k) and h is in the pre- asymptotic range given by 
k^h > 1, then A/12 < kh^ < k^h^ and the -estimate (|5.26|) becomes 

II E - E,,||cG ^Cikh + C2C,tafc'/l'. 

(b) For asymptotic error estimates we refer to 111, section 7.2]. When k^h^ is 
small, it is possible to improve the discrete stability estimates as well as the error 
estimates via the technique of stability- error iterative improvement from fill \2S!\j . 

6. Numerical experiments. Throughout this section, we consider the follow- 
ing Maxwell problem on the unit cube fl = (0, 1) x (0, 1) x (0, 1): 

(6.1) curl curl E - fc^E = in f], 

(6.2) curl E X 1/ - IfcEr = g on T := dn. 

where g is so chosen that the exact solution is E = (e'*^^, e''^^, e'*^^)^. Notice that we 
have chosen X = k for simplicity. 

For any positive integer m, let Ti/m denote the Cartesian mesh that consists of m'^ 
congruent cubes of edge length h = 1/m. We adopt the IPDG method using piecewise 
linear polynomials. We remark that the number of total DOFs of the IPDG method 
on 7i/m is 12m^ which is the about twice of that of the corresponding conforming 
edge element method (EEM) which uses piecewise trilinear polynomials. 

6.1. Stability. Given a Cartesian mesh Th, recall that E/j denotes the IPDG 
solution. Let denotes the trilinear conforming edge element approximation of 

the problem (|6.1I) - (|6.2I) . In this subsection, we use the following penalty parameters 
in the IPDG method ((XTU| : 

(6.3) 70,^ = 70 = 100 and 71,^ = 71 = 0.1 V.F e 
We plot in Figure 16.11 the following two ratios 

III? II IIf,eem|| 

\\'^h\\H{cur\,rh) II h II //(curl, Th) 



l|E||//(curLrfe) l|E||/i-(curI,rh) 

versus A: for /c = 1, 2, • • • , 200 with h = 0.1, 0.05, respectively. It is shown that 

II II //(curl, 7^0 ~ l|E||//(curl,rh) ' 

which is also implied by Theorem 14.21 and Theorem 12.31 The iJ (curl)-norm of the 

edge element solution oscillates for k near 3/h but is still bounded by ||E||^j^j^j.j j-^y 

6.2. Error estimates. In this subsection, we use the same penalty parameters 
as given in (|6.3p . In the left graph of Figure the relative 7J(curl)-error of the 
IPDG solution and the relative iJ(curl)-error of the edge element interpolant are 
displayed in one plot. When the mesh size is decreasing, the relative error of the 
IPDG solution stays around 100% before it is less than 100%, then decays slowly on 
a range increasing with fc, and then decays at a rate greater than —1 in the log- log 
scale but converges as fast as the edge element interpolant (with slope —1) for small 
h. The relative error grows with k along line kh = 1. By contrast, as shown in the 
right of Figure 16.21 the relative error of the finite element solution first stay around 
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200 

(right) 



versus k for k = 1,2, ■ ■ ■ , 200 with h = 0.1, 0.05, respectively. 




l/h 

Fig. 6.2. Left graph: the relative error of the IPDG solution with parameters given in Ij6.3|l 
(solid) and the relative error of the edge element interpolant (dotted) in H{curl)-norm for k = 
5, fe = 10, k = 15, and k = 30, respectively. The dashed line gives reference slope of —1. Right graph: 
corresponding plots for edge element solutions. 

100% but oscillates for large k, then decays at a rate greater than —1 in the log- log 
scale but converges as fast as the edge element interpolant (with slope —1) for small 
h. The relative error of the edge element solution also grows with k along line kh = 1. 

Unlike the error of the edge element interpolant, both the error of the IPDG 
solution and that of the edge element solution are not controlled by the magnitude of 
kh as indicated by the two graphs in Figure l673l It is shown that when h is determined 
according to the "rule of thumb" , the relative error of the IPDG solution keeps less 
than 100% which means that the IPDG solution has some accuracy even for large k, 
while the edge element solution is unusable for large k. We remark that the accuracy 
of the IPDG solution can be further improved by tuning the penalty parameter iji, 
see Subsection 16.31 below. 

Next we verify more precisely the pollution errors. To do so, we recall the def- 
inition of the critical mesh size with respect to a given relative tolerance (cf. [231 
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Fig. 6.3. The relative error of the IPDG solution (left) with parameters given in 116.31 1 and that 
of the edge element solution (right) in H{curl)-norm computed for fc = 2, 4, ■ ■ ■ ,72 with mesh size 
h determined by kh = 2. 

Definition 7.1]). 

Definition 6.1. Given a relative tolerance e and a wave number k, the critical 
mesh size h(k^e) with respect to the relative tolerance e is defined by the maximum 
mesh size such that the relative error of the IPDG solution (or the edge element 
solution) in H(curl)-norm is less than or equal to e. 

It is clear that, if tlie pollution terms are of order fc^/i", then h{k,e) should be 
proportional to k~^/°' for k large enough. Figure which plots h{k,Q.b) versus k 
for the IPDG solution (left) with parameters given in (16.31) and for the edge element 
solution (right), respectively. They all decay at a rate of 0{k~^l'^)^ just like the linear 
FEM for the Helmholtz problem (of. [23]). The results of this subsection indicate that 
both methods satisfy the following pre- asymptotic error bounds (cf. Remark l5.3f a)): 

{\^-^h\\H{ov.rl,TH)^ ||E-E]^™||^^^^^j^^J ^Cikh+C2k^h'^. 

6.3. Reduction of the pollution effect. In this subsection, we show that 
appropriate choice of the penalty parameters can significantly reduce the pollution 
error of the IPDG method. We use the following parameters: 

(6.4) 70,^ = 70 = 100 and i7i,^ = iji = 0.08 + O.Oli VJ" e S^. 

We remark that i7i,j^ is simply chosen from the set {0.01(j3 + qi), —50 ^ p,q ^ 50} to 
minimize the relative error of the IPDG solution in _ff (curl)-norm with 70 = 100 for 
wave number fc = 20 and mesh size h = 1/10. The optimal penalty parameter can 
also be obtained by the dispersion analysis (cf. [1]) and will be considered in a future 
work. 

The relative error of the IPDG solution with parameters given in (|6.4[) and the 
relative error of the edge element interpolant are displayed in the left graph of Fig- 
ure 16.51 The IPDG method with parameters given in (|6.4|) is much better than both 
the IPDG method using parameters given in (|6.3p and the EEM (cf. Figure 16.21 and 
Figure 16. 3p . The relative error does not increase much with the change of k along 
line kh = 1 ior k ^ 30. But this does not mean that the pollution error has been 
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Fig. 6.4. h{k, 0.5) versus k for the IPDG solution (left) with parameters given in 1 16.311 and 
for the edge element solution (right), respectively. The dotted lines give lines of slope —1.5 in the 
log-log scale. 




Fig. 6.5. Left graph: the relative error of the IPDG solution with parameters given in 16. Al 
(solid) and the relative error of the edge element interpolant (dotted) in H{cur\)-norm for k = 
5,k = 10, k = 15, and k = 30, respectively. Right graph: the relative error of the IPDG solution 
with parameters given in 16.41 1 in H{curl)-norm computed for k = 2,4, ■■ ■ ,72 with mesh size h 
determined by kh = 2. 



eliminated. For more detailed observation, the relative error of the IPDG solution 
with parameters given in ()6.4|) in _ff (curl)-norm computed for k = 2,4, •• • ,72 with 
mesh size h determined by fc/i = 2, are plotted in the right graph of Figure 1675] It is 
shown that the pollution error is reduced significantly. 

Figure 1121 plots h{k, 0.5), the critical mesh size with respect to the relative toler- 
ance 50%, versus k for the IPDG method with parameters given in (|6.4I) . We recall 
that h{k, 0.5) is the maximum mesh size such that the relative error of the IPDG so- 
lution in iJ(curl)-norm is less than or equal to 50%. The decreasing rate of h{k, 0.5) 
in the log-log scale is less than —1.5, which means that the pollution effect is reduced. 



For more detailed comparison between the continuous interior penalty finite el- 
ement method (CIP-FEM) and the FEM, we consider the problem (pH - lp?^ with 
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10' ^ , 




10° 10' 10' 

h 

Fig. 6.6. h{k, 0.5) versus k for the IPDG method with parameters given in l|6.4|l . The dotted 
line gives a line of slope —1.5 in the log-log scale. 

wave number k = 36. The real parts of E^3,(0.5, 0.5, z) with parameters given in (I6.4p 
(left, sohd), E™M(Q 5^Q_5^^-) {vigU, sohd), and E^(0.5, 0.5, ^) (dotted) with mesh 
sizes h= 1/18 and 1/36 are plotted in Figure UTT] Here 'Ehx, '^fx^, and E^; are the x 
components of the IPDG solution, the edge element solution, and the exact solution, 
respectively. The shape of the IPDG solution is roughly same as that of the exact 
solution for h = 1/18 and matches very well for h = 1/36. While the edge element 
solution has a wrong shape for h = 1/18 and z > 0.5 and has a correct shape for 
h = 1/36 but suffers an apparent phase error. 

Table 16.11 shows the numbers of total DOFs needed for 50% relative errors in 
i/(curl)-norm for the edge element interpolant, the IPDG solution with parameters 
given in (|6.4p . and the edge element solution, respectively. The IPDG method needs 
less DOFs than the EEM does for k ^ 10 and much less for large wave number k. 



k 


10 


20 


30 


40 


50 


Interpolation 


1,764 


12,168 


33,048 


79,488 


141,288 


IPDG 


2,592 


20,736 


69,984 


187,500 


393,216 


EEM 


2,688 


45,600 


249,900 


876,408 


2,398,488 



Table 6.1 

Numbers of total DOFs needed for 50% relative errors in H{curl)-norm for the edge element 
interpolant, the IPDG solution with parameters given in 1 16.411 . and the edge element solution re- 
spectively. 
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Appendix A. Proof of (|5.19p . The proof follows the same lines as those given 
in [m pages 502-505] and in [IHEl]- Let 

Uh = {vh e H\n); v^,\k e Vif e %} 

be the iJ-'^-conforming linear finite element space. It follows from (|5.ip that E — E/i 
is discrete divergence-free, that is, 

Notice that e V;i n H(curl, il), we have the following discrete Helmholtz decom- 
position of 

(A.l) ^l = Wh + Vrh, 

where rh e Uh and w/i e n H(curl, fl) is also discrete divergence-free. It is easy 
to check that 

l|Vr,i||L2(f2) ||*^||L2(a), W^^hh^in) ^ ||*^||L2(n)- 

Then from |17j, Lemma 7.6] and on noting that the domain 51 is convex, there exists 
w e H^(ri) such that w • n = on F and 

(A. 2) curl w = curl w,,, divw = 0, l|w/i - w||l2(o) < curl *^||L2(n)- 
Thus, it follows from the identity 

(E - E,, ^l)^^ = (E - E,, w,)^^ = (E - E,, w, - w)^^ + (E - E,, w)^ 

that 

(E- E,„*^)„ 

(A.3) ^ \\wh - w||l2(o) + j|w||L2(n) 

||E - Eh||L2(o) 

< h\\ curl*^||L2(j7) + ||w||L2(a). 
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The first term on the right-hand side of (|A.3[) can be bounded as follows: 

(A.4) h\\ curl*^J|L2(n) = h\\ curl(Eft - E + E - E,, + - *^)||L^(n) 

^h\\E,, - E||H(curi,o) +h\\E- EhWoG + h\\ curl(*;, - *^)||l2(7^o < h^TZ{E). 

where we have used (|5.4p . (|5.10l) . and (|5.17l) to derive the last inequality. 

To estimate ||w||l2(j2), we appeal to a duality argument to be described next. Let 
z be the solution of the following auxiliary problem: 

(A. 5) curl curl z + z = w in il, curlz x u = on F := dil. 

Noting that ft is convex, the above problem attains a unique solution z e H^(curl, ft) 
and satisfies the following regularity estimate (cf. |13[ 117) ) 

(A.6) l|zllHi(curi,n) ^ l|w||L2(n)- 

Define sesquilinear forms 

yl(u, v) := (curl u, curl v)o + (u, v)j2 Vu, v e V, 
^/i(u,v) := 6h(u,v) + (u,v)fi Vu,veV. 

Let z^ e V/j n H(curl, il) and z/j e V/j denote the edge finite element approximation 
and the IPDG approximation to z, respectively, that is, 

A{-Vh,z'fJ = (v,i, w)o Vv/i e Y,, n H(curl, fl), 
Ah{vh,Zh) = (v,,,w)o Vv,i e Vh. 

It can be shown that there hold the following estimates (cf. (|5.4I) ): 

l|z - Z^i|H(curl,0) ^ ^ l|z||Hi(curl,n) ^ l|w||L2(n) , 
|||z-Zhll|DG, l||z-Z/,|l|z5G ^ ^(l + 7l)^l|z||Hi(curhn) ^ h {1 + \\mv\\l2 ^a) ■ 

Since 

Mhin) = ^(w, z) = A(w, z - zl) + A(w, z^), 

on noting that w,Wh e H(curl, 57), from (jA.2[) and (|A.4p . we have 

^(w, z - z^) = ^(w - w/i, z - z^) = (w - w/,, z - zl)n 

< h\\ curl*^||L2(o)/il|wl|L2(o) < /^^7^(E)!|w||L2(o). 

On the other hand, 

^(w,z^) = (curl w, curl z^)n + (w,z^)n = (curl curl z^)^^ + (w,z^)o 

= M^K) + (w - (wft + Vrh),zl)^ = A(*^, z^) + (w - w/,, z^)n. 

From the definitions of A, Ah and b/i, we get 
A{^l,zl) = A^m^zl) +\,h{^lzl) 

= Ah{E - Eh,zl) + A(E, - E,z^,) + Ahin - + ^Mn^^l) 

= Ah{B - E,„ z^J + Ahi^l - ^k, zl) + iMn,< - z)- 
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Therefore 

A{w, zl) =Ah{'E - Eh, zl-z) + AhCE - E^, z) + Ah{^l - - z + z - z^) 

+ M^l - ^h, Zh) + i Ji(*^, - z) + (w - w,,, zl)n. 

Since 

Ah{'E-Eh,z) = (E-E,„w)a, - z/0 = (*?,-*h,w)n, 

we have from Lemma 15.11 and the local trace inequality, 

A{w,zl) < |||E-E,j|||z3G|||Zh-z|||DG+ ||E-E,j||L2(n)||w||L2(n) 

+ - *hlllDG(ll|z^ - z|||£,G + ll|z - ZhWlDc) + - */JlL2(n)||wi|L2(n) 
+ 71 II curl*^J|L2(n)(|| curl(z^, - z)||l2(o) + h\\ curl(z^ - z)\\H^rH)) 
+ ||w - wh||L2(a)||z^||L2(n) 
^ I|w||l2(ji) (^h (1 + 7i)^|||E - E,,|||cG + ||E - E,J|l2(ji) + h{l + -fi)^\\\^l - ^hIWdg 

+ -^hh^n) +71^11 curl*^||L2(o) + h\\ curl *^||l2(o)) • 

Moreover, from (|5.17p . (|5.14p . 70 > 1, and the local trace inequality, we get 

111*^ - ^hlWoG ^ (1 + 71)^11 CUrl(*,, - *DllL2(rO + W^hWoG + 11*^ - */i||L2(n) 

< (l + 7i)U7e(E). 
Thus, it follows from (|5.14p . (|A.4|) . ()5.17p . and the above estimate that 

A(w,z^J <(l + 7i)/i27e(E)||w||L2(n). 
Then we obtain the following estimates for ||w||L2(f2-): 

l|w||L2(n) ^ (1 + 7l)/l'7^(E), 
which together with (jA.3[) and (IA.4I) implies that (|5.19p holds. The proof is complete. 



